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ABSTRACT 

We present the results of N-body simulations of planetary systems formation in 
radiatively-inefficient disc models, where positive corotation torques may counter the 
' rapid inward migration of low mass planets driven by Lindblad torques. The aim of 

LlJ , this work is to examine the nature of planetary systems that arise from oligarchic 

growth in such discs. We adapt the commonly-used Mercury-6 symplectic integrator 
by including simple prescriptions for planetary migration (types I and II), planetary 
atmospheres that enhance the probability of planetesimal accretion by protoplanets, 
O ' gas accretion onto forming planetary cores, and gas disc dispersal. We perform a suite 

of simulations for a variety of disc models with power-law surface density and tempera- 
ture profiles, with a focus on models in which unsaturated corotation torques can drive 
outward migration of protoplanets. In some models we account for the quenching of 
corotation torques that arises when planetary orbits become eccentric. Approximately 
half of our simulations lead to the successful formation of gas giant planets with a 
broad range of masses and semimajor axes. Giant planet masses range from being 
approximately equal to that of Saturn, up to approximately twice that of Jupiter. The 
2^ ' semimajor axes of these range from being ~ 0.2 AU, up to ~ 75 AU, with disc models 

that drive stronger outward migration favouring the formation of longer-period giant 
planets. Out of a total of 20 giant planets being formed in our simulation suite, we 
' obtain 3 systems that contain two giants. No super-Earth or Neptune-mass planets 

were present in the final stages of our simulations, in contrast to the large abundance 
of such objects being discovered in observation surveys. This result arises because of 
rapid inward migration suffered by massive planetary cores that form early in the disc 
life time (for which the corotation torques saturate), combined with gas accretion onto 
massive cores that leads them to become gas giants. Numerous low mass planets are 
formed and survive in the simulations, with masses ranging from a few tenths of an 
^ . Earth mass, up to ~ 3 Earth masses. Simulations in which the quenching of corotation 

torques for planets on modestly eccentric orbits was included failed to produce any 
giant planets, apparently because Lindblad torques induce rapid inward migration of 
planetary cores in these systems. 

We conclude that convergent migration induced by corotation torques operating 
during planet formation can enhance the growth rate of planetary cores, but these often 
migrate into the central star because corotation torques saturate. Outward migration 
of planetary cores of modest mass can lead to the formation of gas giant planets at 
large distances from the central star, similar to those observed recently through direct 
imaging surveys. The excitation of planetary eccentricities through planet-planet scat- 
tering during oligarchic growth may quench the effects of corotation torques, however, 
such that inward migration is driven by Lindblad torques. 

Key words: planetary systems, planets and satellites: formation, planet-disc inter- 
actions, protoplanetary discs 



(N 



1 INTRODUCTION 



* E-mail: p.hellary@qmul.ac.uk 



Observations of extrasolar planets are providing a pic- 
ture of a highly diverse population of bodies orbiting 
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main sequence stars. At the extreme ends of the distri- 
bution, there exist very short-period low- mass rocky plan- 
ets such as CoRoT-7 b ijLeeer et all 120091 ) and Kepler- 10b 
(Bata ma et al.l 1201 ll ). and very long-period massive gas 
iant planets detected in recent years by direct imagin g 



Marois et al1l20ld : iKalas et al.ll2008l : iLaerange et~ai1l2010l ) . 
In addition, there have bee n discoveries of short- period 
hot-Jupiters such as 51-Pegb dMayor fc Quelozlll995h . hot- 
Neptunes such as Gliese 436b (iButler et al.l 12004). and 
super- Earths such as Gliese 581c (|Bonfils et al.ll2005l ). Mul- 
tiple planet systems are common, examples being the five 
planet system orbiting the star 5 5 Cancri consisting of gas 
giants and Neptune-mass bodies l|Fischer et al.l l2008l. and 
the recently reported Kepler- 11 s ystem, consisting of six 
nearly coplanar low-mass planets l|Lissauer et al.ll201ll) . A 
question that needs to be addressed is whether or not a par- 
ticular model of planet formation, such as the core-accretion 
model, can explain this broad diversity, appealing to variety 
in the planet forming environment to explain the range of 
observed systems. Or is it the case that quite different phys- 
ical processes are operating on different length and/or time 
scales within protoplanetary discs to form planets with very 
different characteristics. For example, core- accretion operat- 
ing on long-times scales at relatively close locations to the 
central star to form short-period systems, and disc gravita- 
tional fragmentation occurring at large radii on short-time 
scales to form long-period giant planets. 

To begin addressing this question in detail, it is nec- 
essary to construct global models of planetary formation 
that allow for the formation of multiple planet systems 
with a diversity of masses and orbital elements (semima- 
jor axes, eccentricities etc). In this paper, we present the 
results from global models, that are based on the oli- 
garchic growth scenario for planet formation, that have been 
const ructed using a symplectic N-body code (Ch ambers! 
1999), in conjunction with simplified prescriptions for the 
gas disc model, planetary migration, capture of planetes- 
imals, gas-envelope accretion, and disc dispersal on Myr 
time scales. One of our main objectives in this work is 
to examine how our current understanding of migration of 
low-mass planets influences the formation of planetary sys- 
tems, with particular emphasis on the corotation torques 
in radiatively-inefficient discs ( Paardekooper et al.l l20ld: 



iPaardekooper. Baruteau. fe Kiev 20 111 ). As we anticipate 
that this is the first paper in a series that will examine this 
issue, the prescriptions we have adopted in this initial study 
for a number of physical processes, such as gas accretion, are 
necessarily very simplified. They serve the useful purpose, 
however, of enabling N-body simulations to be performed 
of planetary system formation that lead to a diversity of 
outcomes, and these can be used to guide future model de- 
velopments. 

There is a substantial body of previous work that 
has examined the role of migrati on in the formation of 
planet s using N-body simulations. iPapaloizou fc Larwoodl 
(2000) examined planetary growth through planet-planet 
collisions using N-body simulations combined with pre- 
scri ptions for migration and eccentricit y /inc lination damp- 
ing. McNeil. Duncan, fc Levisonl |2005l ) and iDaisaka et. al.l 
(2006) examined the effect of type I migration on ter - 
restrial planet formation, and iFogg fc Nelson! |2007l. 120091 ) 
examined the the influence of type I migration on the 



formation of terrestrial pl anets in the presence of mi - 
grating Jovian-mass planets. iTerquem fc Papaloizoul (|2007l ) 
examined the formation of hot-super-Earths and hot- 
Nept unes using N-body simulatio ns with type I migra- 
tion. iMcNeil fc Nelson! l|2009l . |2010| ) carried out large-scale 
simulations of oligarchic growth to explore the formation 
of systems of hot-Neptunes and super-Earths, such those 
around the stars Gliese 581 and HD698433, using a novel 
symplectic integrator with multiple timesteps. An alter- 
native to these approaches has been p lanetary popula- 
tion s ynthesis modelling, as presented by llda fc Lira (2008. 
l2010l).lMordasini. Alibert. fc Benzl ll2009al)lMordasini et al l 
l|2009bl ). and iMiguel. Guilera. fc Bruninil (|201ll ). These 



monte-carlo approaches have significant advantages in being 
able to cover a very broad range of parameter space, allowing 
meaningful statistical comparisons with observational data 
to be undertaken. The computational efficiency also allows 
complex models of planetary atmospheres and gas accretion 
to be incorporated, as presented by the Mordasini et al work. 
Accurate treatment of planet-planet gravitational interac- 
tions are difficult to include in these models, however, such 
that predictions about planetary system multiplicity, orbital 
eccentricities and inclinations are no t a natural outco me of 
the models (we note that models by llda fc Lin! ((2010) now 
include a simplified treatment of planet-planet interaction 
dynamics). 

This paper is organised as follows. In Sect.[2]we present 
the numerical methods and the physical model. In Sect. [3] 
we present the initial conditions for the simulations. Results 
are described and analysed in Sect.|4j comparisons are made 
to observations in Sect. [5] and a discussion and concluding 
remarks are provided in Sect. [6] 



2 METHOD 

In the following subsections we give details about our phys- 
ical model and the numerical methods we employ. 



2.1 Gas disc model 

To limit the parameter space covered by the simulations, 
we consider only disc models that can provide outward mi- 
gration when unsaturated corotation torques are included. 
The conditions under which outward migration occurs are 
discussed in later sections, but as a rule of thumb we find 
that the temperature radial power-law index, /?, must be ap- 
proximately 0.25 larger than the surface density power-law 
index, a. 

The gas surface density is given by the power-law ex- 
pression 



E 9 (R, t) = E 9 (1 AU) j exp {-t/r disc ) (1) 

where the factor exp (—t/rdisc) mimics the dispersal of the 
gas disc by viscous evolution and photoevaporation on a 
time scale defined by Tdisc- The volume density of gas is 
then 



p(R,Z,t) = 



T,(R,t) 
2^H 



exp(-z 2 /2H 2 ) 



(2) 
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where H is the local disc scale height. The disc temperature 
is also given by a power-law function of radius 



T(R) = T(1AU) 



R 
1AU 



(3) 



A disc with power-law density and temperature profiles also 
has a power-law entropy profile. The associated power-law 
index is given by 



C = (a + /3) - 07 



(4) 



where 7 is the usual ratio of specific heats, here taken to be 
7 = 7/5. The isothermal sound speed is 



m H p 



(5) 



where fee is the Boltzmann constant, mu is the mass of a 
hydrogen atom and p is the mean molecular weight (here 
assumed to equal 2.4). The disc scale height is given by 



H = c s Q}. 



(6) 



where Qk is the Keplerian angular velocity. The angular 
velocity of the gas is given by 



(7) 



2.2 Opacities 

We take the opacity n to be always equal to the Rosseland 
mean opacity, and we take the temperatu re and dens i ty de - 
pendence to be given by the formulae of iBell et all (| 19971 ) 
below 3730 K and bv lBell fc Lml (|l994h above this value: 



«[cm 2 /g] = < 



2 J.-0.01 
rp-1.1 

5 x io 4 r _1 ' 5 

0.1 T 0,7 

2 x 1C 

0.02 T - 8 

2 x 10 81 pT~ 24 

10" 8 p 2/3 T 3 
10 -36 1/3 r io 



2.3 Disc solid component 



T < 132 K 
132 K sC T < 170 K 
170 K sC T < 375 K 
375K sC T < 390K 
390 K sC T < 580 K 
580K sC T < 680K 
680 K < T < 960 K 
960K sC T < 1570K 
1570 K < T < 3730 K 
3730 K < T < 10000 K 



(8) 



The disc solid component is composed initially of protoplan- 
ets and planetesimals (that we model as a computationally 
feasible number of 'superplanetesimals' of much larger mass 
than real planetesimals, but with an assumed radius equal to 
that of realistic planetesimals (10 km) such that they expe- 
rience the appropriate gas drag force.) Protoplanets are ini- 
tially spaced by 10 mutual Hill radii, and planetesimals are 
scattered throughout the disc such that the total solids con- 
tent follows the surface density power-law prescribed for the 



gaseous component. As in lThommes et all (|2003h . planetes- 
imals are distributed according to a Rayleigh distribution 
and have RMS values of the eccentricity e = 0.01 and incli- 
nation i — 0.005 radians, respectively. The surface density 
of solids is enhanced beyond the snow line, whose position 
flsnow is determined by the location where the temperature 
falls below 170 K. The snow line discontinuity is spread over 
a distance ~ 1 AU: 



E s ,o(i?) = <^E 1 + (E 2 -E 1 ) 



1 / R — Rsnow \ 1 \ { R 

2 V 0.5 AU J 2} J V 1AU 

(9) 

The surface density enhancement due to t he snowline 
(E2/E1) = 30/7.1 as in iThommes etH] <|2003h . Planetesi- 
mal densities are set at 3 g/cm 3 throughout the disk. Proto- 
planet densities are set at 3 g/ cm 3 inside the snowlin e and 
1.5 g/cm 3 beyond, as defined bv lThommes et al.1 (2003). The 
mass of the protoplanets at t = is m p = 0.02 Me, and the 
mass of the superplanetesimals is 0.004 M®. 



2.4 Aerodynamic drag 

For kilometre-sized planetesimals, aerodynamic drag pro- 
vides an efficient source of eccentricity and inclination damp- 
ing. We apply gas drag to all bodies in the sim ulation in the 
form of Stokes' drag law (|Adachi et al.lll976h . 

-Fdrag = m p I 3 P D I VlclVrcl (10) 

V 8p P r p J 

where p is the local gas density, p v is the density of the 
planetesimal, r p is the radius of the body and Cd is a di- 
mensionless drag coefficient (here taken to be unity). 



2.5 Capture radii enhancement due to 
atmospheric drag 

If a protoplanet becomes large enough to accumulate a 
gaseous atmosphere, the gas drag acting on planetesimals 
passing through this atmosphere has the effect of increas- 
ing the effective capture radius. We use the p rescription 
described in section 2.5 of llnaba fc Ikomal i|2003ft to model 
this effect (see their Eq. 17 and appendix A). In this model, 
a planetesimal that is within the Hill sphere of the proto- 
planet, and located a distance r c from the centre of the pro- 
toplanet, will be captured if its physical radius is less than 
?"crit given by the expression 



3 v 2 el + 2Gm p /r c p(r c ) 



(11) 



2 v 2 cl + 2Gm p /r H p P 

Here p(r c ) is the local density of the protoplanet atmosphere, 
p p is the density of the planetesimal, th is the protoplanet 's 
Hill radius and « re i is the relative velocity between the two 
bodies. 

The atmosphere model requires us to specify the lumi- 
nosity of the planet. We assume that this is equal to the 
gravitational energy released by incoming planetesimals 

Gmpdrn, 
r p at 

We monitor the accretion rate of solids experienced by 
protoplanets in our simulations, and use this to determine 
the accretion luminosity. In order to smooth out the stochas- 
tic nature of planetesimal accretion, we calculate and use the 
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0.25 



0.15 




0.05 



Figure 1. Effective planetesimal capture radius enhancement due 
to atmospheric drag for 10 km planetesimals and with various 
planet luminosities. Solid lines correspond to L p = 10 — 8 Lq; 
dot-dashed lines correspond to L p = 10 — 5 Lq. 



average luminosity of a protoplanet over temporal windows 
of 200 local orbits (or 4000 years, whichever is smaller). We 
limit the calculated luminosity to within the range 10" 9 to 
10" ' 



The linaba fc Ikomal (|2003h model assumes that the con- 
tribution to the gravitating mass from the atmosphere is 
negligible compared to that of the solid core. In order to 
avoid an overestimation of the capture radius of larger plan- 
ets, we limit the effective capture radius of a planet to be a 
maximum of 1/20 of a planet's Hill radius for these planets. 
The transition is smoothed using the expression 

, - 30M ff 



0.5 - 0.5tanh 



0.5 + 0.5tanh 



5M® 
, - 30M a 

5M ffi 



'"at n 



0.05r H 



(13) 



where rapture is the effective capture radius, r atm os is the 
enhanced capture radius due to the atmosphere and ru is 
the Hill radius. Figure [1] shows the effective capture radius 
as a function of planet mass and luminosity, including the 
above smoothing procedure. 



2.6 Type I migration in radiatively-inefRcient 
discs 



We include the 
our simulations 
mulae presented 



migration of low 
by implementing 



planets in 
torque for- 

by Paardekooper et all (120101 ) and 



mass 
the 



iPaardekooper. Baruteau. fc Klevl ( 20111 ). These formulae 



describe how the various torque contributions vary as 
the planet mass and local conditions in the disc change. 
Specifically, corotation torques depend sensitively on the 
ratio of the horseshoe libration time scale to either the 
viscous or thermal diffusion time scales. 

There are two basic contributions to the corotation 
torque: the vorticity-related corotation torque and the 
entropy-related corotation torque. In an inviscid two di- 
mensional disc, the vortensity (ratio of vorticity to surface 
density) is a conserved quantity on streamlines. Fluid el- 
ements undergoing horseshoe orbits in the presence of a 



planet therefore conserve this quantity. For power-law sur- 
face density profiles with indices greater (less negative) than 
—3/2, there is a negative radial vortensity gradient, and 
the exchange of angular momentum between an embed- 
ded planet and disc material as the fluid follows horse- 
shoe streamlines generating a positive torque on the planet 
l|Goldreich fc Tremainel |l979). In the absence of viscous 
diffusion, material undergoing horseshoe orbits eventually 
phase mixes because of the varying horseshoe orbit time 
scales, erasing the vortensity profile in the corotation region 
and saturating the corotation torque (i.e. switching it off). 
The action of viscous stresses can desaturate the horseshoe 
torque by maintaining the vortensity gradient across the 
horseshoe region, and this occurs optimally when the vis- 
cous diffusion time scale across the width of the horseshoe 
region is approximately equal to half the horseshoe libra- 
tion time. The presence of horseshoe streamlines inevitably 
means that the associated horseshoe torque is a non linear 
effect (because horseshoe orbits are not present in a lin- 
ear theory), usually referred to as horseshoe drag (|Wardl 
Il99ll ; iMassetl l200ll ; IPaardekooper fc Papaloizoul l200St ). As 
the viscosity is increased above its optimal value the vorten- 
sity on streamlines begins to be modified significantly as 
the fluid undergoes horseshoe U-turns. For large enough 
viscosity the vorticity-related corotation torque eventually 
approaches th e smaller value obtained from linear theory 
(|Massetll2002r ). This arises when the viscous diffusion time 
is shorter than the horseshoe U-turn time. 

A similar process occurs for the entropy-related corota- 
tion torque, but in this case the controlling parameter is the 
thermal diffusion time scale instead of the viscous one. Op- 
timal corotation torques are obtained when both the ther- 
mal and viscous diffusion time scales across the width of 
the horseshoe region are equal to approximately half the 
horseshoe libration time. It should be noted that, in addi- 
tion to thermal diffusion, viscosity is required to desaturate 
the entropy-related horseshoe torque. This is because ma- 
terial trapped in the horseshoe region in an inviscid disc 
constitutes a closed system that can only exchange a finite 
quantity of angular momentum with the planet. Viscosity 
is required to couple this region with the rest of the disc, 
such that exchange of angular momentum can desaturate 
the corotation torque. 

For simplicity of im plementation we adopt the ap- 
proxim ation suggested bv lLvra. Paardekooper. fc Mac Low! 
(2010) and assume that the thermal and viscous time scales 
in the disc are equal. For a disc in thermodynamic equilib- 
rium, where the heating is provided by viscous dissipation 
and local cooling is via blackbody radiation, this is a rea- 
sonable assumption to make. 

Based on the above discussion, the torque experienced 
by a low mass planet embedded in a disc depends on the 
Lindblad torques (originating from the excitation of density 
waves at Lindblad resonances), and a weighted sum of the 
vorticity-related horseshoe drag, the entropy-related horse- 
shoe drag, the vorticity-related linear corotation torque, and 
the entropy-related linear corotation torque. These torque 
contributions are given as follows: 
The Lindblad torque is 

T LR = (J^j [-2.5 - 1.7/3 + 0.1a] (14) 
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the vorticity-related horseshoe drag is 



scale 



r I r ° 

1 VHS = I 

7 



l.U--a 



the entropy-related horseshoe drag is 



7.9- 



7 



the vorticity-related linear corotation torque is 



" 7 1 I 



the entropy-related linear corotation torque is 



Tlect = ( — - 

7 



2.2 



— ) c 

7 



(15) 



(16) 



(17) 



(18) 



In these above expression 7 is the ratio of specific heats, and 
To is given by 



2q2 



(19) 



where a p is the planet semimajor axis, and a subscript 
'p' denotes that quantities should be evaluated at the 
orbital location of the planet. In order to obtain the 
correct total torque as a function of the thermal and 
viscous diffusion coefficients we combined the above in- 
dividual torque expressions into the following formula 
l|Paardekooper. Baruteau. fc Klevll201ll ): 

r to t = Tlr + |rvHS-F^Gy + TehsF v FdVGvGd (20) 

+ r LVC T(i - k v ) + r LEC T(i - k u ){\ - k a )}e 

where the functions G v , Gd, F v , Fd, K v and Kd are re- 
lated either to the ratio between the viscous/thermal diffu- 
sion time scales and the horseshoe libration time scale, or to 
the ratio of the viscous/thermal diffusion time scales and the 
horseshoe U-turn time scale. The factor E, that multiplies 
all terms that can contribute to the corotation torque, allows 
for the fact that corotation torques may be strongy attenu- 
ated when the planet has a finite eccentricity, such that it 
undergoes radial excursi ons that are larger th an the width 
of the horseshoe region (|Bitsch fc Kle^koicD . To account 
for this effect we define E according to 

E = (1 -tanh(e/xs). (21) 
where the dimensionless horseshoe width is given by 



1.1 

-vl/4 



(22) 



q — m p /M* and h = H/R. Note that for most simulations 
we set E = 1, but for a subsample of our runs (labelled as 
'E') we use Eq. ([21"]) to define E. 

The horseshoe libration time is given by tub = 
87r/(3f2 p :Es), and the viscous diffusion time scale across the 
horseshoe region is given by t v — (x B a p ) 2 /u, where v is 
the viscous diffusion coefficient. Similarly the thermal dif- 
fusion time scale is given by id = (x B a p ) 2 / D, where D 
is t he thermal diffusion coefficient (d e fined below). Follow- 
ing [^ardekoorjer]^aruteau3k!i£le3 HqUI), we define two 
parameters that determine the relation between the ther- 
mal/diffusion time scales and the horseshoe libration time 



Pv 



16 tv_ 

27 tu b ' 



3 V 2ttv 

which we refer to as the viscous diffusion parameter 



Pd 



2-kD 



3 tub ' 



(23) 



(24) 



which we refer to as the thermal diffusion parameter Note 
that v and D are assumed to be equal in this work, and are 
complicated functions of radial position in the disc because 
of the functional form used to define the opacity in Eq. (JSJ. 
These diffusion parameters are used to define the following 
functions 



(25) 



F d = 



[l + (pv/1.3)*] 
1 

[l + bd/1.3) 2 ]. 



(26) 



Using the viscous diffusion parameter p v we also define the 
following functions 



16 ^457r\3/4p3/2 



16 I 45tt \ < 
25 V 8 / 



if p v < x/8/(457r) 



G v 



K v = 



| - 2 |(4^) 4/3 Pv 8/3 ifpv>y87(45^) 



16 / 457^3/4 3/2 
25 I 28 ) Pv 



ifpv < V8/(45tt) 



(27) 



(28) 



.l-2%(4^) 4/3 Pv 8/3 ifp v >y287(4M 

thermal diffusion 
ctions 

16 (4^)3/4^3/2 ^ < ^ J —- 



Using the thermal diffusion parameter pd we define the fol- 
lowing functions 



Gd 



Kd = 



16 /45tt\3/4 3/2 
25 I 28 I p d 



if p d < v/28/(45tt) 



(29) 



(30) 



2.6.1 Thermal and viscous diffusion 

Radiative diffusion in the disc causes the thermal energy per 
unit volume, e, to evolve according to 

de 



dt 



= -V.F 



(31) 



where the radiative flux in the radial direction (across the 
horseshoe region) may be expressed as 



Fr = ~ 



<ia r c T s dT 



(32) 



3 up dr 

Here a r is the radiation constant and c is the speed of light. 
Noting that e = P/(j — 1) and P = kspT '/ '(fj,mn), and 
assuming that p is locally constant, we obtain the diffusion 
equation governing temperature evolution 



(33) 
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where the diffusion coefficient, D, is given by 
4a r c T 3 (7 - \)jim H 



D 



3 up 2 



kv 



(34) 



We set the viscous diffusion coefficient equai to the thermai 
diffusion coefficient for the purpose of determining the ievei 
of saturation of corotation torques (y = D). 

2.6.2 Eccentricity and inclination damping 

To damp the inclinations of proto planets we us e d the pre- 
scription given in Appendix A of lDaisaka et. all (2006): 



Fidamp ,z — nflp 



M, 



M, 



(35) 

where A c z = -1.088 and A% = -0.871. 

To damp eccentricities we used a simple time scale 
damping formula given by 

0.5 (ve — vk) 



F e damp,r 

where 

tedamp - 



Fedar, 



^edamp 



edamp 



(36) 



Mr. 



-1 / n \ -4 2\ -1 

CLp\ bp \ / 2-igLLp 



fip 1 (37) 



This prescription was adopted rather than using the eccen- 
tricity damping forces prescribed in iDaisaka et. al.1 (120061 ) 
because we found that they could generate significant jitter 
in the acceleration experienced by the planets in disc mod- 
els with strong radial temperature gradients, where H/r be- 
comes small near the disc outer edge. The formula based on 
the time scale argument produced smoother results, appar- 
ently because it is based on an orbit-averaging procedure 
rather than capturing the instantaneous force experienced 
by a planet around its orbit. 

2.7 Gas envelope accretion 

As protoplanets grow through mutual collision and planetes- 
imal accretion they are able to accrete a gaseous envelope 
from the surrounding disc, and may eventually become gas 
giant planets. To model gaseous envelope accretion, we have 
implemented a very approximate scheme by calculating fits 
to the res ults of ID giant planet fo rmation calculations pre- 
sented by iMovshovitz et all (120111 ) . Working in time units 
of Myr and mass units of Earth masses, the gas accretion 
rate is given by 



dm s , 
dt 



5.5 



TgeO 



where we define r gc by the expression 
T gc = 9.665M C "^J C 2 . 
and Tgeo is given by the expression 



TgcO 



e*p(^) 



(38) 



(39) 



(40) 



This procedure allows the planet core to grow due to ac- 
cretion of solids after envelope accretion has commenced, 
and allows the rate of envelope accretion to adjust to the 
changing core mass. It is well known from the studies of 
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Figure 2. Gas accretion onto giant planet cores for 3, 5 and 10 
Mfg cores against time at 5 AU in a disc with no migration or 
planctesimal accretion. 



IPollack et aD (|l996h and others that the rate of gas envelope 
accretion is a sensitive function of the core mass, increasing 
as the core mass increases. Figure [5] shows the gas envelope 
mass evolution in the absence of planetesimal accretion and 
migration for planets with fixed core mass. These are very 
similar to the results of detailed I D giant planet formatio n 
calculations displayed in figure 1 of lMovshovitz et al.l (|201lh . 
Although we have implemented the above equations for gas 
accretion numerically, we note that they have the analytic 
solution 



,(t) = -5.5 In 1 



9.665 



(41) 



Ideally, we would like to incorporate full ID models of giant 
planet formation in our N-body simulations, such that the 
gas envelope accretion can respond to the changing planetes- 
imal accretion rate and changing conditions in the disc, but 
such an approach is computationally prohibitive at present. 
Our simplified approach is highly efficient and provides a 
reasonably good approximation to detailed core nucleated 
accretion models, enabling us to add a vital missing compo- 
nent to our N-body simulations. 

The amount of gas that can accrete rapidly onto a form- 
ing giant planet is constrained by the availability of gas in 
the local feeding zone. We allow giant planets to accrete gas 
using the procedure described above until the envelope mass 
approaches the isolation mass, defined to be the gas mass in 
the feeding zone. W e approximate the feeding zone width to 
be (|Lissauerlll993l ) 

r c = 2\/3r H (42) 

leading to the following expression for the gas isolation mass 
of the planet 



111 . ).,, = / '2:: s '... li dR. 

I a-2V3r H 



(43) 



During the growth of the planet, it can begin to open 
a gap through nonline ar tidal interaction with the disc 
|Lin &: Papaloizoull 19861 ) and for typical disc parameters this 
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occurs around a Jovian mass en et al ]|!999l ; lKlevlll999l ; 
iNelson et all 120001 ). Consequently, if the isolation mass ex- 
ceeds the Jovian mass, we limit the mass of the planet that 
can be obtained during runaway gas accretion to be the Jo- 
vian mass. 

Once the runaway cut-off mass has been reached, the 
gas accretion rate switches to the rate that viscosity can 
supply mass to the planet from the gas disc, 



dt 



(44) 



where v is the local disc viscosity given by 
v = a v H 2 Q{R) 



(45) 

where a v is the viscous parameter (set to 10 -3 for the pur- 
pose of this calculation). Note that this value for the kine- 
matic viscosity is not the same as that obtained by assuming 
the thermal and viscous diffusion coefficients are equal, as is 
done to determine the magnitude of the corotation torques 
acting on a planet (see Sect. 12.6")) . However, the value of a v 
adopted for the purpose of determining the viscously-driven 
mass accretion rate is similar to that used in many previ- 
ous studies of disc-planet int eractions jBrvden et al.lfl999l ; 
iKlevll 19991 ; [Nelson et alJl200C| ). and produces viscous accre- 
tion r ates within the range obse rved to occur onto T Tauri 
stars ( Sic ilia- Aguilar et aUl2004 ). 



2.8 Type II migration 

For massive planets, the migration changes from being of 
type I to being of type II as gap formation occurs. Under 
these circumstances the planet migrates inward on a time 
scale equal to the local viscous evolution time, t„, provided 
that the planet mass is smaller than the local disc mass. 
For more massive planets the migration slows down due 
to the inertia of the planet (and is ultimately determined 
by the time over which the viscous flow in the disc deliv- 
ers a mass of gas co mparable to that of the planet to the 
planet orbital radius (|lvanov. Papaloizou. fc Polnarevll 19991 ; 
ISver fc Clarklll995h ). 

The viscous evolution time is r v = R 2 /(Su), where we 
use Eqn.[45]to calculate u, and we apply the following torque 
in the type II migration regime 



r n = 



m p j p 



1 + 



(46) 



where m p is the planet mass, j p is the specific angular mo- 
mentum, a p is the planet semimajor axis and E p is the disc 
surface density at the planet's semimajor axis location. We 
transition smoothly between the type I and type II migra- 
tion regimes using the following expression 



r cft : = r„s + r I (i-B) 



(47) 



where r e ff is the torque applied during the transition, Ti is 
the type I torque, and the transition function B is given by 



B = 0.5 + 0.5tanh 



65M ff 



15M e 



(48) 



This form for B was adopted to allow the transition to type 
II migration to begin for m p = 30Me , and for the transition 
to be complete for planets with mass m p = 100M® , in broad 
agreement with the results from analytic considerations 



llWardl ll997T ) and numerical simulations (|D'Angelo et al.l 
2003). In the type II regime, eccentricities and inclinations 
are damped on a time scale equal to tv/10. 



3 INITIAL CONDITIONS 

Our simulations were performed u sing the Mercury-6 sym- 
plectic integrator (|Chambers| [l999) . modified to include the 
physics described in Sect. 2. In order to model feasibly mul- 
tiple parameter sets over time scales of 3 Myr, our planetes- 
imal disc consists of super-planetesimals (0.004 Me) with 
effective radius of 10km, representing the averaged orbits 
of a much larger number of real planetesimals. We set the 
mass of our protoplanets to be a factor of 5 larger than the 
planetesimals, giving run times of approximately three cpu 
weeks for each simulation. 

To enable broad coverage of the a and /3 parameter set, 
we limited the number of realisations of initial conditions to 
two runs for each parameter choice, with each member of 
the pair differing only by the random number seed used to 
determine initial positions of the planetesimals. Our initial 
suite of simulations included models with enhancements by 
factors of 3 and 5 abo ve the mass of the Minimum Mass So- 
lar Nebula (MMSN) (|Havashilll98lh (models labelled 'M'), 
but we later augmented these with additional models with 
mass enhancement factors equal to 1 and 2 (models labelled 
'R'). We also examined two models where we implemented a 
reduction in the corotation torques for protoplanets that de- 
velop eccentric orbits (discussed in detail in Sect l2.6p . These 
models are labelled 'E\ Test calculations examining the in- 
fluence of the planetesimal capture radii of protoplanets were 
also performed. 

In order to ensure that the disc mass is locally com- 
parable in models with different surface density profiles, we 
normalise the disc masses so that they all have the same 
mass in the region from 2-15 AU that the enhanced MMSN 
discs would have. This resulted in there being 28 protoplan- 
ets, with ~ 4200 and ~ 2500 planetesimals, for mass en- 
hancement factors of 5 and 3, respectively. We limit our 
selection of the a and j3 parameter space to models for 
which outward migration due to corotation torques is pos- 
sible (the conditions for this can be determined by requir- 
ing the entropy-related and vorticity-related horseshoe drag 
terms in Sect. l2.6l to exceed the Lindblad torques). Our sim- 
ulation parameters are detailed in Table [1] 

We set an inner edge to our simulations at 1 AU, and 
any body that migrates inside this boundary, such that its 
semimajor axis is less than 1 AU, is removed from the sim- 
ulation. Information, however, is stored about each body as 
it crosses this boundary, allowing us to follow the longer 
term trajectories of individual planets to determine their 
final stopping location as the gas disc disperses. This pro- 
cedure is referred to as 'single-body analysis' later in the 
paper. 



4 RESULTS 

In this section we begin by describing common behaviour 
seen in many of the similations. We introduce and discuss 
the concept of zero-migration lines, and their role in creating 
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Figure 3. Contour plots showing regions of outward and inward migration in the mass-semimajor axis plane for runs M05A (a), M16A 
(b), M03B (c) and M07B (d). 



convergent migration within a swarm of growing protoplan- 
ets. We also discuss the coupling between the mass growth 
of protoplanets and their migration, and how rapid accre- 
tion by protoplanets can lead to migration into the central 
star. 

We then describe the detailed evolution of a selection 
of individual runs (four runs in total), followed by a sum- 
mary of results across all of the simulations. This includes 
the results of single-body analyses, where we investigate the 
evolution of bodies lost beyond the inner edge of the sim- 
ulations (these are treated as isolated bodies, and so the 
analyses are limited in their ability to provide accurate pre- 
dictions about the nature of short-period systems). 

Finally, we discuss briefly the effects of protoplanet ec- 
centricity on the collective evolution of the system, and 
present results in which the strength of corotation torques is 
attenuated when a planet's eccentric orbit induces a radial 
excursion that is larger than the horseshoe width. 

Throughout this section, we define a gas giant as being 
a planet that has undergone runaway gas accretion, i.e. the 
sharp increase in mass shown in Fig. [2] This corresponds to 
a mass of approximately 30 M®. 



4.1 Common behaviour 

4-1.1 Migration lines 

Consider a planet orbiting in a protoplanetary disc with 
power-law surface density and temperature profiles. If the 
planet sits in the inner regions of the disc with high sur- 
face density and opacity, such that the horseshoe libration 
time is significantly shorter than the thermal/viscous diffu- 
sion time across the horseshoe region, then the corotation 
torques will saturate and be inoperable. The planet will mi- 
grate inward rapidly as its evolution will be determined by 
Lindblad torques. 

Consider the same planet orbiting much further out 
in the disc where the surface density and opacity are sub- 
stantially reduced, such that the horseshoe libration time is 
much longer than the thermal/ viscous diffusion time. The 
disc-planet system is now close to the locally isothermal 
limit, such that corotation torques will be close to their lin- 
ear value (see Sect. 12. 6] ). The migration will again be inward 
because of the dominance of the Lindblad torques, but at a 
reduced rate because of the contribution of positive corota- 
tion torques. 

There exists an intermediate radial location in the 
disc where the surface density and opacity allow the ther- 
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Table 1. Simulation parameters. 



Simulation 


/cnh 


a 


P 




A'4oiid 


Q 


3110W 


M01A, M01B 


5 


0.5 





75 


173 


1 


95 


M02A, M02B 


3 


0.5 





75 


104 


1 


95 


M03A, M03B 


5 


0.5 


1 




173 


1 


65 


M04A, M04B 


3 


0.5 


1 




104 


1 


65 


M05A, M05B 


5 


0.5 


1 


25 


173 


1 


49 


M06A, M06B 


3 


0.5 


1 


25 


104 


1 


49 


M07A, M07B 


5 


0.5 


1 


5 


173 


1 


39 


M08A, M08B 


3 


0.5 


1 


5 


104 


1 


39 


M09A, M09B 


5 


0.75 


1 




173 


1 


65 


M10A, M10B 


3 


0.75 


1 




101 


1 


65 


M11A, M11B 


5 


0.75 


1 


25 


173 


1 


19 


M12A, M12B 


3 


0.75 


1 


25 


104 


1 


19 


M13A, M13B 


5 


0.75 


1 


5 


174 


1 


39 


M14A, M14B 


3 


0.75 


1 


5 


104 


1 


39 


M15A, M15B 


5 


1 


1 


25 


170 


1 


19 


M16A, M16B 


3 


1 


1 


25 


107 


1 


19 


M17A, M17B 


5 


1 


1 


5 


170 


1 


39 


M18A, M18B 


3 


1 


1 


5 


107 


1 


39 


M19A, M19B 


5 


1.25 


1 


5 


173 


1 


39 


M20A, M20B 


3 


1.25 


1 


5 


101 


1 


39 


R01A, ROIB 


2 


0.5 


1 


25 


69.6 


1 


19 


R02A, R02B 


1 


0.5 


1 


25 


36.6 


1 


49 


R03A, R03B 


2 


0.5 


1 


5 


69.6 


1 


39 


R04A, R04B 


1 


0.5 


1 


5 


36.6 


1 


39 


R05A, R05B 


2 


0.75 


1 


25 


69.6 


1 


19 


R06A, R06B 


1 


0.75 


1 


25 


36.1 


1 


49 


R07A, R07B 


2 


0.75 


1 


5 


69.6 


1 


39 


R08A, R08B 


1 


0.75 


1 


5 


36.1 


1 


39 


EOIA, EOIB 


5 


0.5 


1 


25 


173 


1 


19 


E02A, E02B 


3 


0.5 


1 


25 


104 


1 


49 




Time in years 



Figure 4. Migration lines showing convergent behaviour for 1 
and 2 Me planets in a disc with initial conditions as in M05A. 

mal/viscous diffusion time to be approximately equal to the 
horseshoe libration time. The corotation torque (horseshoe 
drag) will be close to its maximum value here, and will pos- 
sibly drive strong outward migration of the planet if the en- 
tropy gradient in the disc is steep enough. As the planet mi- 
grates outward, however, the local disc surface density and 
opacity decrease, decreasing the thermal diffusion time, and 
reducing the efficacy of the positive corotation torque. Even- 
tually the planet reaches a location where the corotation and 




Figure 5. Migration lines showing convergent behaviour for plan- 
ets with varying mass in a disc with initial conditions as in M03B. 



Lindblad torques exactly cancel, such that the planet stops 
migrating. We refer to this location as the zero-migration 
line, and these are stable positions in the disc for planets to 
reside. 

Given that the horseshoe libration time is shorter for 
more massive planets, the zero-migration lines of heavier 
planets are located further out in the disc where the ther- 
mal diffusion times are shorter. Heavier planets that form in 
inner disc regions will need to migrate out past lower mass 
bodies to reach their zero-migration lines, leading to conver- 
gent migration for protoplanets with different masses. Fur- 
thermore, protoplanets with the same mass try to migrate 
to the same location in the disc. In principle, this should 
increase the rate of collisional planetary growth. 

The behaviour described above is illustrated in the con- 
tour plots shown in Fig. [3] which display the value of the 
total torque in units of To/7 (defined in Eq. I19[). as a func- 
tion of planetary mass and orbital position. The four pan- 
els correspond to the initial models M05A, M16A, M03B 
and M07B that are described in Table. 1. Regions coloured 
red correspond to strong inward migration (migration domi- 
nated by Lindblad torques) . Regions coloured dark blue cor- 
respond to strong outward migration, and lightly coloured 
regions correspond to slow or zero migration. In general, 
rapid outward migration is favoured in discs with relatively 
flat surface density profiles and steep temperature profiles. 

In a steady-state disc, a planet of fixed mass that mi- 
grates to its zero-migration line should stay there. Disc dis- 
persal, however, leads to a locally reducing surface density 
and opacity, progressively shifting the zero- migration line in- 
ward. Consequently, as the disc disperses, a planet sitting on 
a zero-migration line drifts inward on the gas disc dispersal 
time scale. This behaviour is shown in Figs. [4] and [5] which 
show the migration trajectories of planets of different mass 
in the two disc models M05A and M03B that are dispersing 
on time scales of 1 Myr (similar migration traje ctories are 
shown in lLvra. Paardekooper. fc Mac Low! (|201Cf O. Planets 
of a given mass starting at different locations tend to migrate 
outward and eventually join the same migration line, which 
then moves inward as the disc disperses. The behaviour of 
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Figure 6. Contour plots showing regions of outward and inward 
t=2,000,000 years (c) and t=3,000,000 years (d). 

the contours shown in the top left panel of Fig. [3] under the 
action of disc dispersal are shown in Fig. [6] It is clear that 
as the disc becomes increasingly optically thin, only heavy 
planets can sustain outward migration (unless they become 
too massive and enter the type II migration limit because of 
gap formation). 



4-1.2 Influence of mass growth 

A planet undergoing mass growth while sitting on a zero- 
migration line should migrate outward as it accretes, since 
zero-migration lines for heavier planets lie at larger distances 
from the central star. However, if the mass growth rate of 
the planet exceeds the outward migration rate due to coro- 
tation torques, then a very different fate awaits it. Consider 
a planet of mass ~ 0.5 Me sitting at ~ 2 AU in the top left 
panel of Fig. [6] Rapid growth of the planet up to 10 in 
less than 1 Myr will put the planet in the regime of rapid in- 
ward migration, as its trajectory in the figure will be almost 
vertical, moving it out of the blue region and into the red 
one. Very rapid growth of planets, therefore, may not lead 
to strong outward migration but instead may cause planets 
to migrate rapidly into the central star. The timing of the 
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migration for run M05A at t=0 years(a), t=l,000,000 years (b), 



growth of planets is crucial in determining their long-term 
survival. 

4.2 Individual runs 

4.2.1 Run M05A 

Run M05A has an initial disc mass equivalent to 5 times the 
MMSN, a = 0.5 and /3 = 1.25. The magnitude and sign of 
migration torques (t = 0) are shown in the top left panel of 
Fig. and the effects of disc dispersal are demonstrated in 
Fig- El It is clear that planets with masses in the range 0.2 ^ 
m p ^ 1 M(f), initially located in the disc region 1-10 AU, 
can undergo strong outward migration. Growth of planets 
to masses of a few may lead to outward migration over 
distances of tens of AU. 

The time evolution of planet masses (top panel), semi- 
major axes (middle panel) and eccentricities (bottom panel) 
are shown in Fig. [7] During the first 0.3 Myr, we observe that 
three planets grow in mass rapidly, and undergo outward mi- 
gration to ~ 10 AU. The mass growth occurs as a result of 
planetesimal accretion and planet-planet inelastic collisions, 
and the rapid growth is assisted by convergent migration 
within the protoplanet swarm and by the gas envelopes that 
form within the planet Hill spheres. When the planet masses 
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Figure 7. Evolution of the masses, semimajor axes and eccentricites of all protoplanets simulation M05A. 



exceed ~ 20 Mq , however, their migration direction changes 
and they undergo very rapid inward migration through the 
planetary swarm and interior to 1 AU, the inner boundary 
of the simulation. During the rapid inward migration, there 
is very little accretion by these bodies, but they temporarily 
excite the eccentricities of the other bodies in the system 
(see the bottom panel of Fig. [7] between 0.2 - 0.3 Myr). 

Between the times 0.3 - 0.5 Myr, we observe that three 
bodies grow to masses larger than 1 Mj . The largest of these 
grows to ~ 4 Mq and migrates outward rapidly to ~ 30 AU, 
the location of its zero-migration line. We refer to this as 
"planet A". A second planet ("planet B") grows to a mass 
~ 3 Mq by 0.5 Myr, and migrates out to its zero-migration 
line at ~ 20 AU. A third planet ("planet C") reaches a mass 
of ~ 2 M® at 0.5 Myr and migrates out to 10 AU. Although 
the protoplanet/planetesimal disc of solids is truncated at 
15 AU in the initial simulation set-up, outward migration 
of planets and gravitational scattering transports planetesi- 
mals into the outer disc where they are accreted by planets 
A and B (planet C continues to reside within the original 
planetesimal disc). Gas accretion ensues once these bodies 
exceed 3 M®, and we see their masses grow smoothly up to 
20 - 30 M® between the times 0.5 - 2 Myr. During this time 
the zero-migration lines move inward (observe the modest 
inward migration in the middle panel of Fig. [7] between 0.5 
- 2 Myr), but continued mass growth helps to counterbal- 
ance this effect, and prevents substantial inward migration. 
At approximately 2 Myr, planets A and B undergo rapid 
gas accretion and grow to become Jovian-mass giant plan- 
ets (further gas accretion occurs at the viscous-supply rate). 
The rapid mass growth induces dynamical instability be- 
tween planets A and B, causing them to undergo a period 



of gravitational scattering and eccentricty growth (bottom 
panel of Fig. [7]) . The scattering eventually causes planets B 
and C to collide at 2.1 Myr (when planet C has a mass of 24 
M®), leaving two giant planets on eccentric, non-crossing or- 
bits with semimajor axes of 12 and 55 AU. The inner planet 
has a total mass of 374 M®, and a solid core mass of 27.6 
M®, at the end of the simulation. The outer planet has a 
total mass of 352 M®, and a solid core mass of 10.1 M®. 

During the formation of the outer gas giant planets, 
between time 0.5 - 2.5 Myr, only modest planetary growth 
occurs in the inner syst em. An inner resonant convoy, simi - 
lar to those discussed bylMcNeil. Duncan, fc Levi son (2005) 
and lCresswell fc Nelson! <|2006t ) migrates interior to 1 AU by 
~ 1.6 Myr, driven by a more rapidly migrating 0.5 M® body. 
This leaves behind two planets that grow to become ~ 3 and 
0.4 M® before migrating interior to 1 AU at ~ 2.5 Myr. 



4.2.2 Run M16A 

Run M16A has a disc mass equivalent to 3 times the MMSN, 
a — 1 and ft — 1.25. The migration behaviour at t — is il- 
lustrated by the contours displayed in the top right panel of 
Fig. [3 and it is clear that outward migration is considerably 
weaker in this model than in the previously described run 
M05A. Furthermore, the steepness of the outward migra- 
tion region as one moves to higher planet masses indicates 
that the radial extent of outward migration is also reduced 
relative to model M05A. Placed in the initial disc model, a 
planet with m p < 1 M® orbiting at 1 AU will not undergo 
outward migration at all, but will instead migrate inward 
only. Rapid planetary growth is therefore expected to result 
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Figure 8. Evolution of the masses, semimajor axies and eccentricities of all protoplanets in simulation M16A. 

in much of the solid disc material migrating in toward the 
star. 

The evolution of the planetary masses (upper panel), 
semimajor axes (middle panel) and eccentricities (bottom 
panel) are shown in Fig. [8] Protoplanets located initially 
beyond ~ 2 AU in this disc with masses ~ 0.02 M® (the 
inital masses of protoplanets in the model) are expected to 
migrate inward, and looking at the middle panel of Fig. [5] 
we see obvious evidence for this migration occuring within 
the first 0.3 Myr. Looking at the inner part of the system 
during the first 0.5 Myr, we observe two examples of reso- 
nant, inward migrating convoys being established. The first 
to form consists of the six innermost protoplanets in the 
system. All masses of these planets are < 1 M®, except for 
the outermost body, whose mass has grown to ~ 1 M®. 
The more rapid migration of this body drives the inward 
migration of the whole convoy. At a time of ~ 0.4 Myr, we 
see that the next three nearest protoplanets to the central 
star begin to undergo rapid inward migration, and this is 
driven by the formation of a ~ 5 M® body who's progenitor 
protoplanet was located at ~ 4 AU. The growth of this pro- 
toplanet induces rapid inward migration, with the system of 
inner planets being swept interior to 1 AU at t = 0.55 Myr. 

Three planets initially located at ~ 5 AU become phys- 
ically detached from the rest of the system after ~ 0.5 Myr, 
as shown in the middle panel of Fig. [8] These bodies have 
all grown to masses between 2-5 M® within this time. 
The outermost ~ 2 M® body becomes isolated from plan- 
etesimals in the disc such that its mass does not grow after 
0.5 Myr. This isolation occurs in large part because the two 
more massive neighbouring planets accrete the nearby plan- 
etesimals. Having achieved masses in excess of 3 M®, these 



two planets are able to accrete gas. As they do so, they 
sit on their zero-migration lines and undergo slow inward 
migration, where the speed of migration is attenuated by 
the continuing gas accretion and mass growth (the planets 
try to migrate outward to the zero-migration lines for more 
massive planets as they grow, at the same time as the zero- 
migration lines move inward as the gas disc evolves). After 
~ 2.6 Myr, the innermost planet reaches a mass of ~ 30 
M® and undergoes rapid gas accretion to become a Saturn- 
mass gas giant. The rapid mass growth dynamically disturbs 
the system, as observed in the middle and bottom panels of 
Fig. [8] causing the outer 2 M® planet to collide with the 
middle planet. Shortly after, this merged planet undergoes 
rapid gas accretion to also become a Saturn-mass gas gi- 
ant. Saturnian rather Jovian masses are achieved because 
accretion occurs late in the disc lifetime, such that the gas 
isolation mass limits the envelope mass to ~ 100 M®. 

At the end of the simulation we have an inner planet of 
total mass 115 M®, with solid core mass 11 M®, orbiting at 
2.3 AU, and an outer planet with total mass 127 M®, solid 
core mass 8.8 M®, orbiting at 3.1 AU. 

4.2.3 Run M03B 

Run M03B has a disc mass equivalent to 5 times that of 
the MMSN, a = 0.5 and /3 = 1. The migration behaviour 
of this model at t — is illustrated by the lower left panel 
in Fig. O showing that this disc is intermediate between the 
two previous models discussed (M05A and M16A) in terms 
of the strength of outward migration. 

The evolution of the planetary masses (upper panel), 
semimajor axes (middle panel) and eccentricities (bottom 
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Figure 9. Evolution of the masses, semimajor axes and eccentricities of all protoplanets in simulation M03B. 



panel) are shown in Fig. [9] As expected from comparing 
the migration behaviour of the runs M16A and M03B in 
Fig. [3] the initial stages of run M03B show some similarities 
to run M16A. Protoplanets initially located in the region of 
the protoplanet/planetesimal disc between 10 - 15 AU mi- 
grate inward to the region centred around 2-3 AU. The 
inner planets, however, do not show a strong tendency to 
migrate inward (differing in this regard from run M16A), 
and the convergent migration stimulates substantial growth 
within the protoplanet swarm, as seen in the upper panel 
of Fig. [9l where the planet masses are seen to increase dur- 
ing the first 0.5 Myr, and in the middle panel where it is 
clear that collisional growth reduces the number of planets. 
The strong planetary growth, however, also leads to rapid 
inward migration. Bodies that reach masses in excess of 20 
undergo rapid inward migration through the disc of pro- 
toplanets/planetesimals and interior to 1 AU, exciting the 
eccentricities of the remaining planets as they do so. The 
bodies that rapidly migrate through the inner boundary at 
1 AU would hit the star if their long-term evolution were 
followed. 

One of the quickly migrating planets (shown by the up- 
per green line in the top panel of Fig. [9j grows to be massive 
enough (approximately 30 Mj of solids) to undergo rapid 
gas accretion during the inward migration. It reaches its gas 
isolation mass at a mass equal to 114 M®, and transitions to 
type II migration, drifting interior to 1 AU shortly after 0.8 
Myr has elapsed. At this point in time, the planet mass is 
~ 250 M®, and it is undergoing gas accretion from the disc 
at the viscous-supply rate. We have followed the evolution 
of this body after it has traversed the inner boundary of the 
simulation, treating it as an isolated body and neglecting 
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Figure 10. Evolution of mass and semimajor axis of single body 
extension run for the hot Jupiter in simulation M03B starting at 
0.821 Myr. 

its interaction with other bodies in the system (we refer to 
this as single-body analysis). The evolution is displayed in 
Fig. 1101 and we see that the planet reaches a semimajor axis 
of 0.25 AU and has a mass of 524 M® after 3 Myr, making 
it an excellent candidate for a hot Jupiter. 

No other planets grow substantially during the evolu- 
tion of this run. Fig. [9] shows that only a single planet with 
mass ~ 0.35 M® survives beyond 1 AU, coming to rest at a 
semimajor axis of ~ 1.9 AU. 

4.2.4 Run M07B 

The final run we describe in detail is M07B, which has a 
mass equivalent to 5 times the MMSN, a = 0.5 and /3 = 1.5. 
The steep temperature gradient and relatively shallow sur- 
face density gradient allow this disc model to support strong 
outward migration over a large radial extent, as illustrated 
by the contours in the bottom right panel of Fig. [3] This 
plot demonstrates clearly that sub-Earth mass bodies orbit- 
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ing in the vicinity of 1 AU will experience strong outward 
migration, possibily out beyond 100 AU. 

The evolution of planetary masses (top panel), semi- 
major axes (middle panel) and eccentricities (bottom panel) 
are shown in Fig. 1111 The initial growth and outward mi- 
gration of protoplanets in the inner region of the swarm 
leads to strongly convergent migration, and rapid growth of 
a number of bodies up to ~ 10 M® within the first 0.3 Myr. 
As observed in the previously described runs, this leads to 
rapid inward migration of these planets because the horse- 
shoe libration time becomes much shorter than the ther- 
mal/viscous diffusion time for these bodies. After 0.5 Myr, 
there are only three bodies left in the simulation: two plan- 
ets orbiting at ~ 15 AU each with masses ~ 2 M®; one 
protoplanet with mass ~ 0.3 M® orbiting at ~ 2 AU. The 
two outer bodies collide shortly after 0.5 Myr, and the re- 
sulting planet undergoes slow gas accretion, migrating out- 
ward as it does so. After 2.8 Myr it undergoes rapid gas 
accretion and becomes a Jovian mass (319 M®) gas giant 
planet, with a 5.2 M® solid core, orbiting at ~ 33 AU. As in 
model M05A, we find that gas giant planets can be formed 
at large radius from the central star by the migration and 
gas accretion onto a solid core. In both of these models (and 
others not discussed in detail), the mode of formation is one 
in which an initial generation of massive protoplanets form 
and migrate in toward the central star, followed by a sec- 
ond generation of more isolated lower mass cores that can 
migrate out slowly and accrete gas at the same time. Such 
a model may provide a natural explanation for the massive 
planets orbiti ng at large distance from their hos t stars, such 
as HR 8799 |Marois et al.ll201Cf) . Fomalhaut l|Kalas et all 



2008), and Beta Pictoris (|Lagrange et aljfeoicj ). that are 
being discovered by direct imaging surveys. 

4.3 Summary of all runs 

We ran 40 simulations with / ea h = 5 or 3, surface density 
power-law indices in the range 0.5 a ^ 1.25, and tem- 
perature power-law indices satisfying 0.75 ^ /3 ^ 1.5. In 
total 19 gas giant planets were formed in these runs, and 
their properties are summarised in Fig. 1 121 and Table [2] The 
giant planet masses range from 115 to 670 M®, and have 
solid core masses in the range 3.6 - 39 M®. Looking at the 
upper panel of Fig. 1121 we see that most giant planets are 
grouped within the mass range 320 to 400 M®, and this is 
very likely to be an artifact of our gas accretion procedure 
that limits the planet mass obtained during rapid gas accre- 
tion to be the Jovian mass. A more sophisticated procedure 
would be sensitive to local conditions in the disc, and re- 
sult in a broader spread of planet masses, and this is clearly 
one future improvement that we will need to implement in 
our modelling procedure. Nonetheless, we do also obtain gi- 
ant planets outside of this mass range. Run M03B formed 
a 523.8 M® planet, as discussed in Sect. 14.231 due to a gas 
giant forming within the first 0.5 Myr, and subsequently ac- 
creting viscously and migrating via type II migration to its 
final location at 0.25 AU. The heaviest planet formed during 
run M11A, and this was the result of two gas giant planets 
colliding, having each formed at between 20 and 30 AU from 
the central star due to their cores migrating outward. Em- 
ploying a pure hit-and-stick prescription for planetary colli- 
sions, however, probably leads to an overestimate of the final 
mass of this planet. Three planets were formed with masses 
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below the imposed Jovian-mass limit. Run M12B produced 
a 296 planet orbiting at 9.8 AU, and as described in 
Sect. 14.2.21 run M16A produced a pair of Saturn-mass ob- 
jects orbiting at 1.15 and 1.55 AU. These planets formed 
late in the disc lifetime, such that their gas isolation masses 
were below the Jovian mass. 

Most of the giant planets formed at semimajor axes 
substantially beyond 10 AU. Indeed, only 4 out of 19 giant 
planets formed interior to 10 AU. The reason for this is that 
the most common mode of gas giant planet formation in the 
simulations was the formation of a core of modest mass in 
the interior disc, that then migrates outward over large dis- 
tances before accreting a gas envelope. Many massive cores 
formed during the early stages of the disc life times in the 
simulations, and were able to undergo gas accretion. Their 
rapid inward migration, however, prevented them from sur- 
viving. Giant planets that form closer to the star and survive 
tend to be in disc models (M03B and M16A) that generate 
weaker corotation torques. 

There are three simulations that lead to the formation 
of surviving multiple giant planets. Run M05A produces a 
pair of Jovian mass planets orbiting at 11.4 and 53.9 AU, 
as described in Sect. 14.2.11 and M05B also produced a pair 
of Jovian mass objects orbiting at 23.3 and 68.8 AU. Note 
that these runs were identical apart from the random num- 
ber seed used to generate initial conditions. Run M16A pro- 
duced a pair of Saturn-mass planets orbiting at 1.15 and 1.55 
AU. One consequence of this is that almost all giant plan- 
ets formed in the simulations are on circular, non-inclined 
orbits. The only planets with significant eccentricities are 
those in run M05A, where gravitational scattering during 
the formation caused the growth of eccentricity. 

One surprising result to come out of the simulations is 
the lack of correlation between initial disc mass and the fre- 
quency of giant planet formation. Discs with / cn h = 5 formed 
9 giants and those with / cn h = 3 formed 10. This led us to 
question whether less massive discs might be able to form 
gas giants. To examine this, we performed additional simu- 
lations (labelled 'R' in Table [TJ based on the most succesful 
models in the f en h=3 and 5 runs. All barring one failed to 
produce any gas giants. Run R07B, a / e „h=2 disc, did pro- 
duce a single Jovian mass gas giant (shown in Fig. I12p . 

In addition to the giant planets discussed above, the 
simulations also resulted in the formation and survival of 
lower mass bodies beyond 1 AU in the disc. These are shown 
in Fig. 1131 The rapid growth of cores, followed by rapid in- 
ward migration, has the effect of clearing much of the solid 
material from beyond 1 AU, and the outward migration of 
modest sized cores that evolve into gas giants also clears this 
region. Nonetheless, terrestrial mass bodies do form and sur- 
vive in the simulations, although Fig. [T3] shows that these 
tend to be in the lower mass discs. One noticeable and inter- 
esting observation about the simulation results is the lack 
of super-Earth and Neptune mass planets. The rapid for- 
mation of massive cores that undergo fast inward migration 
is a major cause of this (driven by efficient capture of plan- 
etesimals and convergent migration), but a contributing fac- 
tor is the fact that planets with masses greater than 3 Mj 
can begin to undergo gas accretion in our models. A higher 
threshold for gas accretion would probably allow some of the 
giant planets that formed to have maintained lower masses. 
These observations provide a useful guide to the types of 



modifications that the modelling procedure requires in or- 
der to form planets with the same characteristics as those 
which are contained in the extrasolar planet observational 
database. 

4-3.1 Single-body analysis 

The inner edge of our simulations was set at 1 AU. We ran 
single body runs for each object that was lost beyond this 
inner edge so as to identify any bodies that would have be- 
come short period gas giants, and to obtain an estimate of 
the distribution of smaller bodies in this inner region. These 
runs are effectively continuation runs, but with a smaller 
time step size, and a smaller inner boundary at 0.1 AU. It 
is important to note that these single body runs did not 
include the influence of any other planets in the system, 
and did not account for any material that would have been 
present between 0.1 - 1 AU during the early evolution of the 
system. As such, the results merely provide a guide to the 
planets that can survive within this radial range. 

Figure [14] shows a summary of all non-giant planets left 
remaining in the 0.1 to 1 AU region from all the f en h=3 and 
5 models. Objects with masses less than 1 M® are clearly 
more common than those with larger masses because of their 
reduced migration rates. Also, smaller semimajor axes seem 
the more likely outcome. All objects included in this figure 
have masses below 3 and have not been able to undergo 
gas accretion. The only gas giant to survive in the region 0.1 
- 1 AU is the one described already in Sect 14.2.31 

A large number of bodies are lost beyond 0.1 AU, rang- 
ing from the smallest protoplanets all the way up to Jupiter 
sized gas giants, potentially providing the central star with 
a significant enrichment of heavy elements. These bodies are 
summarised in Fig. 1151 which shows the masses of the plan- 
ets as they cross the boundary at 1 AU in the lower panel, 
and their masses as they cross the boundary at 0.1 AU in 
the upper panel. The horizontal axes show the time of loss 
through the boundary at 0.1 AU. It is clear that a number of 
the massive cores that migrate through the 1 AU boundary 
accrete gas and become gas giants, although their masses 
normally reach values between 100 - 200 because the 
gas isolation mass is below the Jovian-mass in the inner disc. 
Type II migration drives them through the boundary at 0.1 
AU. It is also clear that a number of bodies migrate inside 1 
AU with masses in the range 4-10 and grow through 
gas accretion to masses between 30 - 50 M®. Type I migra- 
tion, however, forces these bodies to migrate into the star 
before they can become giants. Their corotation torques are 
saturated, and so rapid inward migration is driven by Lind- 
blad torques. 

Some of the bodies passing through the 0.1 AU bound- 
ary at late times could have survived. We ran extended single 
body runs for the five planets with masses greater than 25 

lost beyond 0.1 AU in the last 500,000 years of simula- 
tion time. Two collided with the central body at 2.75 Myr, 
but the other three survived at 0.086, 0.0428 and 0.016 AU 
with masses 344, 164 and 550 M®, respectively. We have not 
included these results along with the other gas giants since 
their simulation conditions were overly simplified compared 
to the rest. All three bodies entered the 1 AU zone with 
just a few Earth masses, and so would in reality have in- 
teracted with other bodies and planetesimals formed there 
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Figure 12. Summary of total masses, core masses, eccentricity and inclination against semimajor axis for all gas giant planets formed. 



3 
2.5 

2 

°> . _ 
E 

1 

0.5 



2 * Minimum Mass 

3 * Minimum Mass 
5 * Minimum Mass 



1 



4 5 
semimajor axis (AU) 



Figure 13. Summary of masses against semimajor axis for all small surviving planets outside 1 AU. 



which were not modelled. Also, the body ending up at 0.016 make reasonable predictions about the nature of surviving 
AU would most likely have been accreted by the central star short-period planets, 
a short while later. 



The survivability of these giant planets formed in single 
body analyses between 1 and 0.1 AU depends on exactly 
how the gas disc dissipates. The exponential dissipation of 
gas provides a reasonable approximation for the bulk of the 
gas dissipation when it is dominated by viscous evolution 
(Fogg fc Nelson! l2007t ). but at later times the structure of 
a viscously evolving disc that is being photoevaporated by 
UV radiation from the central star changes substantially 
ijClarke. Gendrin. fc Sotomavorjl200ll ) with a low density in- 
ner cavity being formed. Clearly such a model needs to be 
incorporated into the simulation procedure outlined here to 



4.4 Eccentricity modulation of corotation torques 

Recent numerical simulations |Bitsch fc Klevll2010r i indicate 
that corotation torques are substantially reduced in their 
effectiveness when a planet develops an eccentric orbit. In 
particular, we expect the corotation torque to be effectively 
quenched when the radial excursion associated with the ec- 
centric orbit exceeds the width of the horseshoe region. In 
order to simulate this effect, we have run a few simulations 
(labelled 'E' in Table [T]) where the eccentricity modulation 
function in Eq. [21] is switched on. The effect of this was 
as one might expect: growth was significantly stunted com- 
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Table 2. Summary of gas giants formed. 



Simulation 


/enh 


a 


P 




a (AU) 


e 


i (degrees) 


mass (Mgj) 


core mass (Mgj) 


M03B 


5 


0.5 


1 




0.24818 


0.000001 





523.79 


39.39 


M05A 


5 


0.5 


1 


25 


11.39435 


0.125762 


2.7366 


374.36 


27.57 


M05A 


5 


0.5 


1 


25 


53.91049 


0.191585 


1.5949 


352.2 


10.11 


M05B 


5 


0.5 


1 


25 


23.26593 


0.00314 


0.0299 


351.47 


9.51 


M05B 


5 


0.5 


1 


25 


68.79704 


0.011044 


0.0195 


433.61 


19.18 


M06A 


3 


0.5 


1 


25 


54.99131 


0.000684 


0.0026 


369.38 


12.11 


M06B 


3 


0.5 


1 


25 


74.69739 


0.000847 


0.0012 


392.5 


16.1 


M07A 


5 


0.5 


1 


5 


55.91612 


0.000897 


0.0006 


319.86 


3.6 


M07B 


5 


0.5 


1 


5 


32.59897 


0.00098 


0.0057 


319.24 


5.19 


M08A 


3 


0.5 


1 


5 


13.41661 


0.001873 


0.0106 


374.61 


26.13 


M08B 


3 


0.5 


1 


5 


60.55699 


0.000785 


0.0021 


333.38 


13.35 


M11A 


5 


0.75 


1 


25 


24.93238 


0.063072 


0.0053 


669.88 


19.7 


M11B 


5 


0.75 


1 


25 


22.70169 


0.000959 


0.0069 


361.13 


7.79 


M12B 


3 


0.75 


1 


25 


9.87598 


0.000726 


0.0022 


296.43 


6.47 


M14A 


3 


0.75 


1 


5 


27.22201 


0.000179 


0.0046 


323.17 


5 


M14B 


3 


0.75 


1 


5 


49.00695 


0.000497 


0.0036 


328.24 


12.07 


M16A 


3 


1 


1 


25 


1.55415 


0.008245 


0.0008 


114.91 


10.9511 


M16A 


3 


1 


1 


25 


1.15495 


0.007372 


0.0008 


126.85 


8.8343 


M18B 


3 


1 


1 


5 


23.61206 


0.000516 


0.0111 


325.91 


7.95 


R07B 


2 


0.75 
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5 


54.41991 


0.000599 


0.0016 


367.3 
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Figure 14. Summary of masses against semimajor axis for all small surviving planets interior to 1 AU. Note that these data were 
obtained using the single-body analysis described in the text. 



pared to the corresponding runs without this reduction fac- 
tor (runs M05A/B compare to eccentricity damping runs 
E01A/B and M06A/B compare to E02A/B). Nearly all pro- 
toplanets were lost beyond the inner edge by approximately 
1 Myr and most protoplanets were lost within half this time. 
Only one planet survived to run completion out of all four 
of the simulations and its final position is shown in Fig. 1181 
which shows a summary of surviving planets from the 'E' 
runs, as well as those discussed below in which the enhanced 
planetesimal capture radii are switched off. A plot showing 
the time evolution of this particular simulation is given in 
Fig. 1161 where we have plotted the eccentricity in the bot- 
tom panel as e/x 3 . It is clear that planet-planet interactions 
maintain values of e/x s ^ 1 throughout the simulation, un- 
til it is depleted of planets through their inward migration. 
This result suggests that closer investigation of the role of 
planetary eccentricity in regulating the strength of corota- 
tion torques needs to be undertaken, since the modest evi- 
dence we have accumulated suggests that mutual encounters 
between planets may remove the benefits provided by coro- 
tation torques in enhancing the formation and survival of 



planets. Similar conclusions have been reached in a recent 
study by McNeil & Nelson (In preparation) that examines 
the formation of hot Neptunes and super-Earth planets in 
radiatively inefficient discs. 



4.5 Capture radii enhancement switched off 

A common outcome within our simulations has been the 
rapid formation and growth of planetary cores, and their 
subsequent rapid migration inward. One reason for this 
rapid growth is our adoption of an enhanced capture ra- 
dius for planetesimals arising because of a putative gaseous 
envelope settling onto protoplanets during their formation. 
We re-ran the simulations described in Sect. 14.21 without 
enhanced capture radii. Growth was notably slower as ex- 
pected in all four runs. Two runs, however, (corresponding 
to the M03B and M07B non-atmosphere runs) did manage 
to form one gas giant planet each. 

The time evolution of run M03B-NA (non-atmosphere) 
is shown in Fig. 1171 A planet grows slowly to just over 3 Me 
by approximately 500,000 years. It sits in an area of the disc 
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Figure 16. Evolution of masses, semimajor axes and eccentricities of all protoplanets in simulation E02B. 



Planet formation in optically thick discs 19 




0.25 - 
0.2 - 

a, 0.15- 




0.5 1 1.5 2 2.5 3 

Time (years) 10 s 



Figure 17. Evolution of masses, semimajor axes and eccentricities of all protoplancts in simulation M03B-NA. 
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Figure 18. Summary of masses against semimajor axis for all small surviving planets outside 1 AU for both runs where eccentricity 
damping was turned off and where enhanced capture radii due to atmospheric drag were turned off. 
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largely cleared of solid material by other protoplanets and 
slowly accretes gas before eventually undergoing runaway 
gas accretion at 2.5 Myrs and ends up at 2 AU with a mass 
of 126 M®. 



Run M07B-NA forms a gas giant by means of three 1- 
1.5 M® bodies migrating out to large semimajor axes, and 
then merging to form a 3.5 Mj body at 40 AU which slowly 
accretes gas until runaway gas accretion occurs at 2.8 Myrs. 
The planet ends up at 50 AU with a mass of 319 M®. The 
lower mass planets that survived in these runs are shown in 
Fig. 1181 Interestingly, one of these is a Neptune-sized planet. 



5 COMPARISON WITH OBSERVATIONS 

The work presented in this paper is not intended to be 
a serious attempt at planetary popul ation syn t hesis mod- 
elling, unlike th e work presented by llda fc Lml (|200&t ) and 
iMordasini et al.l (|2009bh . Instead it is aimed at exploring 
the consequences of having strong corotation torques oper- 
ating during the oligarchic growth stage of planetary sys- 
tems formation, and understanding how planetary growth, 
migration and planet-planet interactions combine to form 
planetary systems. Nonetheless, it is of interest to examine 
how well the simple models that we have presented compare 
with the observational data on extrasolar planets. 

Figure [19] shows a mass-period diagram with our re- 
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Figure 19. Mass vs semimajor axis plot comparing observed 
cxtrasolar planets with our simulation results. 




1(T 2 10"' 10° 10 1 10 2 10 3 10 4 



semimajor axis (AU) 

suits overlaid on all current observed exoplanets (sourced 
from www.exoplanet.eu). Our shorter period giant planets 
lie well in the range of already detected exoplanets both in 
terms of mass and semimajor axis. Our longer period plan- 
ets, however, lie in an area that is sparsely populated with 
observational results. It is worth noting, however, that this 
region of parameter space is much more problematic for the 
detection of planets, as observations rely on direct imaging 
methods rather than radial velocity or transit observations. 

A clear failing of our results is in the formation of super- 
Earths and Neptune-mass bodies. One reason for this ap- 
pears to be the gas accretion routine that we have adopted, 
that allows planets to accrete gas once their masses exceed 
3 M$. An additional issue is the adoption of atmosphere- 
induced enhanced capture radii for planetesimal accretion 
onto protoplanets. The very rapid growth of planets due to 
this often causes them to migrate rapidly toward the central 
star, an outcome that is reduced in models where enhanced 
capture radii are not included. These are issues that we will 
address in a future publication. 



6 DISCUSSION AND CONCLUSIONS 

Although the models presented in this paper include a broad 
range of physical processes relevant to planet formation 
(migration; planetary growth through mutual protoplanet 
collisions; accretion of planetesimals and gas; planet-planet 
gravitational interactions), we have adopted a number of 
assumptions and simplifications that inevitably affect the 
realism of the simulations and their results. These include: 

(i) Simulation domain. Even though we have modelled a 
relatively wide semimajor axis domain with our initial solid 
matter disc compared to previous N-body work on plane- 
tary formation, it is clear that accurate modelling of discs 
in which significant corotation torques arise requires as wide 
a domain as possible. Protoplanets move significantly in the 
disc with some forming at 2-3 AU and migrating out to 70-80 
AU. Similarly, planets migrate inward and ought to traverse 
the terrestrial planet region which we have not modelled. 
Planets forming in the terrestrial region might also migrate 
out into the regions that we have investigated. In short, the 



migration behaviour observed in the simulations presented 
in this paper indicates that all regions of the disc are cou- 
pled during planet formation, and it is no longer sensible to 
think in terms of a "terrestrial planet region" or a "giant 
planet region". As such, a suitable model would involve a 
domain ranging from as far in as 0.1 AU out to approxi- 
mately 50 AU. Such a simulation is beyond current mod- 
elling techniques because of the required numbers of proto- 
planets and plane t esima l s, eve n for the method presented by 
iMcNeil fc Nelson! l|2009l . |201Ch which utilises multiple time 
steps in a parallel symplectic integrator. Instead, such global 
models of planetary formation are probably going to require 
efficient use of modern GPU technology. 

(ii) Gas disc model. We currently model the gas disc as 
having fixed power-laws in surface density and temperature, 
with the disc mass undergoing exponential decay to mimic 
its viscous and photoevaporative evolution. In reality, the 
disc is heated by the central star and through local vis- 
cous dissipation, and it cools through radiative emission. A 
significant improvement to the model that we are in the 
process of implementing will be to evolve the disc surface 
density and temperature explicity using a 1+1D nume rical 
solution, as described in IPapaloizou fc Terqueml (|l999h . for 
example. This a pproa ch will be similar to that described in 
iFogg fc Nelson] (2009) and references therein, and will al- 
low gap formation and type II migration to be simulated 
directly, along with UV photoevaporation of the disc during 
its final stages of dispersal. 

(iii) Planetary atmosphere model and enhanced capture 
radii. As described in the preceding sections, rapid plane- 
tary growth is assisted by the enhanced accretion of plan- 
etesimals through implementation of a model for planetary 
atmo spheres that increase s the effective accretion cross sec- 
tion (|lnaba fc Ikomall2003h . Although this model works well 
when accretion is dominated by planetesimals, it is probably 
inaccurate when accretion includes giant impacts between 
protoplanets. In particular, a planetary atmosphere can be 
completely liberated from a planet when it is impacted by a 
body whose mass is similar to that of the atmosphere, and 
our implementation of the atmosphere model does not ac- 
count for this effect. The atmosphere model would clearly 
be improved in its accuracy if it responded to giant impacts 
as well as planetesimal accretion rates. 

(iv) Gas envelope accretion. Our model for gas accretion 
during the formation of gas giants is very rudimentary, al- 
though it serves the purpose of allowing gas giant planets to 
form in our simulations. While a ful l accretion model for each 
planet similar to those p resented in iPollack et ail (|l996T ) and 
iMovshovitz et all (|201ll ) would be ideal, this is computation- 
ally beyond the reaches of an N-body code that can model 
planetary systems formation over Myr time scales. However, 
there are improvements that can be made that will allow lo- 
cal conditions in the disc to influence the gas accretion rate 
onto a planet. Coupled with a more sophisticated disc model 
that allows explicit modelling of gap formation and gas ac- 
cretion, such an approach would alleviate the requiement to 
set an artifical upper limit for the planet mass that can form 
through runaway gas accretion. 

We have presented the results of simulations that in- 
clude recent torque formulae for type I migration (including 
Lindblad and corotation torques), with gas envelope accre- 
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tion, enhanced capture radii due to gas atmospheres, and 
planet-planet gravitational dynamics included. We have sur- 
veyed a range of disc models which all allow for outward mi- 
gration driven by corotation torques. The main results that 
we have obtained may be summarised as follows: 

• Convergent migration of protoplanets, and the rapid ac- 
cretion of planetesimals, can cause the rapid growth of plan- 
etary cores to masses in excess of 10 within 0.5 Myr in 
most disc models. This leads to rapid inward migration of 
these bodies, driven by Lindblad torques, when the horse- 
shoe libration time scale becomes significantly shorter than 
the thermal/viscous diffusion time scale and the corotation 
torques saturate. 

• Formation of planetary cores with a few Earth masses 
^ 0.5 Myr after the simulations are initiated can lead to their 
migration into the outer regions of the disc (30 - 50 AU). 
Steady mass growth through gas accretion onto the planet 
can counterbalance the slow inward migration that occurs 
as the gas disc mass reduces, and long-term survival in the 
outer disc can lead to gas giant planet formation there when 
runaway gas accretion ensues. This mode of giant planet for- 
mation was found to be a common outcome in our simula- 
tions, especially those with disc models that generate strong 
outward migration, leading to numerous gas giant planets 
being formed between semimajor axes 10 - 60 AU. Mod- 
els such as these could potentially explain the long-period 
giant planets discovered recently through direct imagin g 
jMarois et ai]|2010l ; iKalas et al.ll2008l ; lLagrange et~alll2010l ) . 

• Out of 40 simulations that used disc models whose 
masses were either 3 or 5 times more massive than the 
MMSN, 19 gas giant planets were formed. Most of these 
are similar in mass to Jupiter (in part because of the gas 
accretion prescription that was adopted in the models), and 
are formed at large distances from the star. Short period 
Jovian mass planets were also formed, however, along with 
a pair of Saturn-mass bodies at intermediate (~ 1 AU) or- 
bit distances. These latter systems were formed in discs that 
generate weaker corotation torques than those that tend to 
generate the longer-period giant planets. 

• Multiplanet systems containing more than one giant 
planet were found to be a rare outcome (3 out of 40 simu- 
lations produced two giant planets each), and this has the 
additional effect of producing systems with very small ec- 
centricities and inclinations due to the low rate of occurence 
of planet-planet scattering events. In fact, the only planets 
to be formed with significant eccentricities were a pair of 
closely separated Jovian-mass objects that underwent sig- 
nificant dynamical interaction. 

• Our simulations completely fail to produce super-Earth 
or Neptune-mass planets. This appears to arise because of 
very rapid inward migration of planets that grow early in the 
disc lifetime and undergo rapid inward migration, combined 
with the switching-on of gas accretion that converts planets 
of intermediate mass into gas giants at later times. Modifi- 
cation of the planetary atmosphere and gas accretion pre- 
scriptions will probably result in more surviving planets with 
intermediate masses. Numerous planets in the Earth mass 
range were formed in the simulations, however. The 'desert' 
of super-Earths and Neptunes is similar to tha t reported 
in the planetary population synthesis models of llda fc Linl 



(2008), and occurs for much the same reasons as theirs (rapid 
gas accretion and migration). 

• Simulations performed where the corotation torque is 
attenuated when planet eccentricities grow to become larger 
than the dimensionless horseshoe width appear to produce 
results quite different from those in which this effect is ne- 
glected (i.e. all the runs described above). In particular the 
growth and survival of planets is reduced because mutual en- 
counters between protoplanets maintains the typical eccen- 
tricities above the critical value for which corotation torques 
diminish. In these latter simulations, no gas giant planets 
were formed at all. Further work is required to establish the 
influence of corotation torques on planet formation via oli- 
garchic growth, where planet-planet interactions maintain 
finite values of the eccentricity. 

Our models demonstrate that strong corotation torques 
can substantially alter the qualitative outcomes of planet 
formation simulations. Even the simplest model of planet 
formation that involves non linear gravitational interaction 
between protoplanets and planetesimals during planetary 
accumulation is by its nature a chaotic process. Given a 
well-defined distribution of initial protoplanet/planetesimal 
masses and orbital elements from which individual forma- 
tion models are drawn, however, an ensemble of such mod- 
els should give rise to a distribution of outcomes with well- 
defined statistical properties. Allowing the set of initial con- 
ditions to be drawn from a range of disc models whose life- 
times and radial profiles of density and temperature also 
have well-defined distributions serves only to modify the dis- 
tribution of outcomes, as does including additional physical 
processes such as type I migration. A corollary of this ar- 
gument is that increasing the complexity of migration pro- 
cesses, as we have done in this paper, also serves to modify 
the distribution of outcomes in a quantifiable manner. Coro- 
tation torques, however, increase the dependency of forma- 
tion outcomes on the details of the underlying disc model 
and microphysical processes such as those that control the 
opacity of disc material. The dependency of migration on 
opacity, turbulent viscosity and other disc properties, and 
the role of migration in shaping planetary system architec- 
tures, increases the need for more refined observations of 
protoplanetary disc properties and improved disc models to 
allow planetary formation calculations to be compared with 
data on extrasolar planetary systems in a meaningful way. 
To summarise: planetary formation is a chaotic process, but 
is deterministic in a statistical sense. Corotation torques do 
not change the validity of this statement, however, their de- 
pendence on detailed disc properties increases the difficulty 
of constructing realistic planetary formation models for com- 
parison with observational data. 

This is the first in a series of papers to examine the 
oligarchic growth of planets under the influence of type I 
migration, including corotation torques. Models that include 
a more sophisticated treatment of the issues raised in Sect. [6] 
will be presented in a forthcoming publication. 
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